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1. Introduction 


The Fermi Gamma-ray Space Telescope, which was launched by NASA in June 2008, is a 
powerful space observatory which studies the high-energy gamma-ray sky Atwood (2009). 
Fermi’s main instrument, the Large Area Telescope (LAT), detects photons in an energy range 
between 20 MeV to greater than 300 GeV. The LAT is much more sensitive than its predecessor, 
the EGRET telescope on the Compton Gamma Ray Observatory, and is expected to find 
several thousand gamma-ray point sources, which is an order of magnitude more than its 
predecessor EGRET Hartman et al. (1999). 


Even with its relatively large acceptance (~2 m? sr), the number of photons detected by the 
LAT outside the Galactic plane and away from intense sources is relatively low and the sky 
overall has a diffuse glow from cosmic-ray interactions with interstellar gas and low-energy 
photons that makes a background against which point sources need to be detected. In 
addition, the per-photon angular resolution of the LAT is relatively poor and strongly energy 
dependent, ranging from more than 10° at 20 MeV to ~0.1° above 100 GeV. Consequently, 
the spherical photon count images obtained by Fermi are degraded by the fluctuations on 
the number of detected photons. This kind of noise is strongly signal dependent : on 
the brightest parts of the image like the galactic plane or the brightest sources, we have 
a lot of photons per pixel, so the photon noise is low. Outside the galactic plane, the 
number of photons per pixel is low, which means that the photon noise is high. Such a 
signal-dependent noise can’t be accurately modeled by a Gaussian distribution. The basic 
photon-imaging model assumes that the number of detected photons at each pixel location is 
Poisson distributed. 


More specifically, the image is considered as a realization of an inhomogeneous Poisson 
process. This statistical noise makes the source detection more difficult, consequently it is 
highly desirable to have an efficient denoising method for spherical Poisson data. 
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Several techniques have been proposed in the literature to estimate Poisson intensity in 
2D. A major class of methods adopt a multiscale bayesian framework specifically tailored 
for Poisson data Nowak & Kolaczyk (2000), independently initiated by Timmerman & 
Nowak (1999) and Kolaczyk (1999). Lefkimmiaits et al. (2009) proposed an improved 
bayesian framework for analyzing Poisson processes, based on a multiscale representation 
of the Poisson process in which the ratios of the underlying Poisson intensities in adjacent 
scales are modeled as mixtures of conjugate parametric distributions. Another approach 
includes preprocessing the count data by a variance stabilizing transform (VST) such as 
the Anscombe Anscombe (1948) and the Fisz Fisz (1955) transforms, applied respectively 
in the spatial Donoho (1993) or in the wavelet domain FryZlewicz & Nason (2004). The 
transform reforms the data so that the noise approximately becomes Gaussian with a constant 
variance. Standard techniques for independant identically distributed Gaussian noise are 
then used for denoising. Zhang et al. (2008) proposed a powerful method called Multi-Scale 
Variance Stabilizing Tranform (MS-VST). It consists in combining a VST with a multiscale 
transform (wavelets, ridgelets or curvelets), yielding asymptotically normally distributed 
coefficients with known variances. The interest of using a multiscale method is to exploit 
the sparsity properties of the data: the data is transformed into a domain in which it is 
sparse, and, as the noise is not sparse in any transform domain, it is easy to separate it 
from the signal. When the noise is Gaussian of known variance, it is easy to remove it 
with a high thresholding in the wavelet domain. The choice of the multiscale transform 
depends on the morphology of the data. Wavelets represent more efficiently regular structures 
and isotropic singularities, whereas ridgelets are designed to represent global lines in an 
image, and curvelets represent efficiently curvilinear contours. Significant coefficients are 
then detected with binary hypothesis testing, and the final estimate is reconstructed with an 
iterative scheme. In Starck et al. (2009), it was shown that sources can be detected in 3D LAT 
data (2D+time or 2D+energy) using a specific 3D extension of the MS-VST. 


To denoise Fermi maps, we need a method for Poisson intensity estimation on spherical 
data. It is possible to decompose the spherical data into several 2D projections, denoise 
each projection and reconstitute the denoised spherical data, but the projection induces some 
caveats like visual artifacts on the borders or deformation of the sources. 


In the scope of the Fermi mission, two of the main scientific objectives are in a sense 
complementary: 


e Detection of point sources to build the catalog of gamma ray sources, 
e Study of the Milky Way diffuse background. 


The first objective implies the extraction of the Galactic diffuse background. Consequently, we 
want a method to suppress Poisson noise while extracting a model of the diffuse background. 
The second objective implies the suppression of the point sources: we want to apply a binary 
mask on the data (equal to 0 on point sources, and to 1 everywhere else) and to denoise the 
data while interpolating the missing part. Both objectives are linked: a better knowledge 
of the Milky Way diffuse background enables us to improve our background model, which 
leads to a better source detection, while the detected sources are masked to study the diffuse 
background. 


The aim of this chapter is to present a multi-scale representation for spherical data 
with Poisson noise called Multi-Scale Variance Stabilizing Transform on the Sphere 
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(MS-VSTS) J.Schmitt et al. (2010), combining the MS-VST Zhang et al. (2008) with various 
multi-scale transforms on the sphere (wavelets and curvelets) Abrial et al. (2007); Abrial 
et al. (2008); Starck et al. (2006). Section 1.2 presents some multi-scale transforms on the 
sphere. Section 1.3 introduces a new multi-scale representation for data with Poisson noise 
called MS-VSTS. Section 1.4 applies this representation to Poisson noise removal on Fermi 
data. Section 1.5 presents applications to missing data interpolation and source extraction. 
Section 1.6 extends the method to multichannel data. 


All experiments were performed on HEALPix maps with nside = 128 Gorski et al. (2005), 
which corresponds to a good pixelisation choice for the GLAST/FERMI resolution. 


2. Wavelets and curvelets on the sphere 


New multi-scale transforms on the sphere were developed by Starck et al. (2006). These 
transforms can be inverted and are easy to compute with the HEALPix pixellisation, and 
were used for denoising, deconvolution, morphological component analysis and inpainting 
applications Abrial et al. (2007). In this chapter, here we use the Isotropic Undecimated 
Wavelet Transform (IUWT) and the Curvelet Transform. 


2.1 The HEALPix pixellisation for spherical data 


Different kinds of pixellization scheme exist for data on the sphere. For Fermi data, we 
use the HEALPix representation (Hierarchical Equal Area isoLatitude Pixellization of a 
sphere) Gorski et al. (2005), a curvilinear hierarchical partition of the sphere into quadrilateral 
pixels of exactly equal area but with varying shape. The base resolution divides the sphere 
into 12 quadrilateral faces of equal area placed on three rings around the poles and equator. 
Each face is subsequently divided into nside2 pixels following a quadrilateral multiscale tree 
structure. (Fig. 1) The pixel centers are located on iso-latitude rings, and pixels from the same 
ring are equispaced in azimuth. This is critical for computational speed of all operations 
involving the evaluation of spherical harmonic transforms, including standard numerical 
analysis operations such as convolution, power spectrum estimation, etc. HEALPix is a 
standard pixelization scheme in astronomy. 


2.2 Isotropic Undecimated Wavelet Transform on the sphere 


The Isotropic Undecimated Wavelet Transform on the sphere (IUWT) is a wavelet transform 
on the sphere based on the spherical harmonics transform and with a very simple 
reconstruction algorithm. At scale j, we denote a;(0, p) the scale coefficients, and d;(0, p) 
the wavelet coefficients, with 0 denoting the longitude and ¢ the latitude. Given a scale 
coefficient a jr the smooth coefficient a j+1 is obtained by a convolution with a low pass filter h j: 
dj+1 = 4j * hj. The wavelet coefficients are defined by the difference between two consecutive 
resolutions : dj,1 = 4; — 4j41. A straightforward reconstruction is then given by: 
J 
ag (0, p) =a; (0, p) +) d;(8, p) (1) 
j=l 


Since this transform is redundant, the procedure for reconstructing an image from its 
coefficients is not unique and this can be profitably used to impose additional constraints 
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Fig. 1. The HEALPix sampling grid for four different resolutions. 


on the synthesis functions (e.g. smoothness, positivity). A reconstruction algorithm based on 
a variety of filter banks is described in Starck et al. (2006). Figure 2 shows the result of the 
IUWT on WMAP data (Cosmic Microwave Background). 


2.3 Curvelet transform on the sphere 


The curvelet transform enables the directional analysis of an image in different scales. The 
data undergo an Isotropic Undecimated Wavelet Transform on the sphere. Each scale j is 
then decomposed into smoothly overlapping blocks of side-length B; in such a way that 
the overlap between two vertically adjacent blocks is a rectangular array of size Bj x B;/2, 
using the HEALPix pixellisation. Finally, the ridgelet transform Candes & Donoho (1999) is 
applied on each individual block. The method is best for the detection of anisotropic structures 
and smooth curves and edges of different lengths. The principle of the curvelet transform is 
schematized on Figure 3. More details can be found in Starck et al. (2006). 


2.4 Application to Gaussian denoising on the sphere 


Multiscale transforms on the sphere have been used successfully for Gaussian denoising via 
non-linear filtering or thresholding methods. Hard thresholding, for instance, consists of 
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Fig. 2. WMAP data and its wavelet transform on the sphere using five resolution levels (4 
wavelet scales and the coarse scale). The sum of these five maps reproduces exactly the 
original data (top left). Top: original data and the first wavelet scale. Middle: the second and 
third wavelet scales. Bottom: the fourth wavelet scale and the last smoothed array. 


setting all insignificant coefficients (i.e. coefficients with an absolute value below a given 
threshold) to zero. In practice, we need to estimate the noise standard deviation 0; in each 
band j and a coefficient wj is significant if |w;| > xo;, where x is a parameter typically chosen 
between 3 and 5. Denoting Y the noisy data and HT, the thresholding operator, the filtered 
data X are obtained by: 

X = @HT,(®'Y), (2) 


where ®! is the multiscale transform ([UWT or curvelet) and ® is the multiscale 
reconstruction. A is a vector which has the size of the number of bands in the used multiscale 
transform. The thresholding operation thresholds all coefficients in band j with the threshold 
Aj = 0; 

J Ji 
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= Curvelets on the Sphere 


Rudgelet Transform on the Sphere (RTS) 


Fig. 3. Principle of curvelets transform on the sphere. 
3. Multi-scale transforms on the sphere and Poisson noise 
3.1 Principle of the Multi-Scale Variance Stabilizing Transform on the Sphere (MS-VSTS) 


In this section, we propose a multi-scale representation designed for data with Poisson 
noise. The idea is to combine the spherical multi-scale transforms with a variance stabilizing 
transform (VST), in order to have a multi-scale representation of the data where the noise 
on multi-scale coefficients behaves like Gaussian noise of known variance. With this 
representation, it is easy to denoise the data using standard Gaussian denoising methods. 


VST of a Poisson process 


Given Poisson data Y := (Y;);, each sample Y; ~ P(A;) has a variance Var[Y;] = A;. Thus, the 
variance of Y is signal-dependant. The aim of a VST T is to stabilize the data such that each 
coefficient of T(Y) has an (asymptotically) constant variance, say 1, irrespective of the value 
of A;. In addition, for the VST used in this study, T(Y) is asymptotically normally distributed. 
Thus, the VST-transformed data are asymptotically stationary and Gaussian. 


The Anscombe Anscombe (1948) transform is a widely used VST which has a simple 


square-root form 
T(Y) :=2VY+3/8. (3) 
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We can show that T(Y) is asymptotically normal as the intensity increases. 
D 
T(Y) — 2VA——>N (0,1) (4) 
A — +00 


It can be shown that the Anscombe VST requires a high underlying intensity to well stabilize 
the data (typically for A > 10) Zhang et al. (2008). 


VST of a filtered Poisson process 


Let Zj := i h[i]Yj—; be the filtered process obtained by convolving (Y;); with a discrete filter 


h. We will use Z to denote any of the Z;’s. Let us define % := y;(h[i])* for k = 1,2,-++. In 
addition, we adopt a local homogeneity assumption stating that Àj- i = A for all i within the 
support of h. 


We define the square-root transform T as follows: 
T(Z) := b- sign(Z + c)|Z + c|", (5) 


where b is a normalizing factor. It is proven in Zhang et al. (2008) that T is a VST for a 
filtered Poisson process (with a nonzero-mean filter) in that T(Y) is asymptotically normally 
distributed with a stabilized variance as A becomes large. 


The Multi-Scale Variance Stabilizing Transform on the Sphere (MS-VSTS) consists in 
combining the square-root VST with a spherical multi-scale transform (wavelets, curvelets...). 
3.2 Wavelets and Poisson noise 


This subsection describes the MS-VSTS + IUWT, which is a combination of a square-root VST 
with the IUWT. The recursive scheme is: 


IUWT aj = hj—1 *aj_-4 
dj = Aaj} =j 
(6) 
MS-VSTS a; h;_ * AJ 
=$ { j j—-1*#j—1 : 
+IUWT | dj = Tj-1(4j-1) — T;(a;) 


In (6), the filtering on a;_; can be rewritten as a filtering on ag := Y, i.e., a; = hi) x ag, where 
8 j 8 j 

AY) = hj_y * +++ * hy * hg for j > land hn) = ô, where 6 is the Dirac pulse (ô = 1 ona single 

pixel and 0 everywhere else). Tj is the VST operator at scale j: 


T;(aj) = bY sign(a; +.) \/|a; + c0]. (7) 


Let us define ol! Vie L; (hU) [i])*. In Zhang et al. (2008), it has been shown that, to have an 
optimal convergence rate for the VST, the constant c) associated to hU) should be set to: 
G) G) 
Fc) PE ir en a (8) 
8r ) an! ) 
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The MS-VSTS+IUWT procedure is directly invertible as we have: 


J 
a9(8,9) = Ty! [nen + Eaj (8, p). (9) 
j=l 


Setting bÜ) := sign(tl! )) / ee ) |, if A is constant within the support of the filter. hO), then we 
have Zhang et al. (2008): 


A => +00 ari (10) 
ol) (hG, i) 
n Oe 


where (.,.) denotes inner product. 


This means that the detail coefficients issued from locally homogeneous parts of the signal 
follow asymptotically a central normal distribution with an intensity-independant variance 
which relies solely on the filter h and the current scale for a given filter A. Let us define T) the 


stabilized variance at scale j for a locally homogeneous part of the signal: 


=T H P r 

TE Z (Av), nO) 11 

eme e G0) ` an 
4t; AT) 2T T 


To compute the Tj) bO) cG dll ) we only have to know the filters hl). We compute these 
filters thanks to the formula a; = hO) x ag, by applying the IUWT to a Dirac pulse ag = ô. 


Then, the h are the scaling coefficients of the IUWT. The Tj) have been precomputed for a 
6-scaled IUWT (Table 1). 


Wavelet scale j 
1 


0.484704 
0.0552595 
0.0236458 


2 
3 
4 0.0114056 
5 0.00567026 


Table 1. Precomputed values of the variances gj of the wavelet coefficients. 


We have simulated Poisson images of different constant intensities À, computed the IUWT 
with MS-VSTS on each image and observed the variation of the normalized value of 0; j) 
((œŒ(j))simulated / (o j) theoretical) as a function of À for each scale j (Fig. 4). We see that the 
wavelet coefficients are stabilized when A 2 0.1 except for the first wavelet scale, which 
is largely noise. In Fig. 5, we compare the result of MS-VSTS with Anscombe + wavelet 
shrinkage, on sources of varying intensities. We see that MS-VSTS works well on sources of 


very low intensities, whereas Anscombe does not work when the intensity is too low. 


www.intechopen.com 


Poisson Noise Removal in Spherical Multichannel Images: Application to Fermi Data 169 


sigma(1) sigma(2) 
12 r~ r T 12 T T T 
1.0 1.0 
+ 
0.BF 4 0.8F 4 
+ 
0.6F 4 os 4 
o4ab 4 oal 4 
o2 E 4 0.26 4 
0.0L 1 1 1 fi 4 0.0L 1 1 1 1 
o 2 4 6 8 10 o 2 4 6 8 10 
lambda lambda 
sigma(3) sigma(4) 
12 — T T Lap T T T 
Lop | tok + % 
ob 4 o8 F + 
o.s E 4 o.s E 4 
O4b 4 oal 4 
o2 E 4 o2 E 4 
0.9 1 1 1 0.0 , 1 1 , 
o 2 8 10 o 2 4 6 8 10 
lambda lambda 
sigma(5) 

12 — T T 

106 + 

ope ed 

o.s E 4 

oF 4 

o.2F 4 

v.o L L 1 1 1 

o 2 4 6 8 10 


lambda 


Fig. 4. Normalized value ((ce (j)) simulated / (7 )) theoretical) of the stabilized variances at each 
scale j as a function of À. 


3.3 Curvelets and Poisson noise 


As the first step of the algorithm is an IUWT, we can stabilize each resolution level as in 
Equation 6. We then apply the local ridgelet transform on each stabilized wavelet band. 


It is not as straightforward as with the IUWT to derive the asymptotic noise variance in the 
stabilized curvelet domain. In our experiments, we derived them using simulated Poisson 
data of stationary intensity level A. After having checked that the standard deviation in 
the curvelet bands becomes stabilized as the intensity level increases (which means that the 
stabilization is working properly), we stored the standard deviation 0; ; for each wavelet scale 
j and each ridgelet band / (Table 2). 
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Fig. 5. Comparison of MS-VSTS with Anscombe + wavelet shrinkage on a single face of the 
first scale of the HEALPix pixelization (angular extent: 7t /3sr). Top Left : Sources of varying 
intensity. Top Right : Sources of varying intensity with Poisson noise. Bottom Left : Poisson 
sources of varying intensity reconstructed with Anscombe + wavelet shrinkage. Bottom Right 
: Poisson sources of varying intensity reconstructed with MS-VSTS. 


0.348175 
2| 0.230621 | 0.248233 | 0.196981 


3| 0.0548140 |0.0989918} 0.219056 
4| 0.0212912 |0.0417454 |0.0875663 | 0.20375 
5|0.00989616 |0.0158273 |0.0352021 |0.163248 


Table 2. Asymptotic values of the variances g; ; of the curvelet coefficients. 
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4. Application to Poisson denoising on the sphere 
4.1 MS-VSTS + IUWT 


Under the hypothesis of homogeneous Poisson intensity, the stabilized wavelet coefficients d j 
behave like centered Gaussian variables of standard deviation Tj). We can detect significant 
coefficients with binary hypothesis testing as in Gaussian denoising. 


Under the null hypothesis Ho of homogeneous Poisson intensity, the distribution of the 
stabilized wavelet coefficient dj|k] at scale j and location index k can be written as: 


1 
Pld) = Fare PAK); (12) 
The rejection of the hypothesis Ho depends on the double-sided p-value: 
iS io exp(—x2/202)dx (13) 
4 V270; J|d;lk]| P ae 


Consequently, to accept or reject Ho, we compare each |dj|k]| with a critical threshold xo;, 
K = 3,4 or 5 corresponding respectively to significance levels. This amounts to deciding that: 


e if |dj[k]| > xo;, d;[k] is significant. 
e if |dj[k]| < xo;, d;[k] is not significant. 


Then we have to invert the MS-VSTS scheme to reconstruct the estimate. However, although 
the direct inversion is possible (Eq. (??)), it can not guarantee a positive intensity estimate, 
while the Poisson intensity is always nonnegative. A positivity projection can be applied, but 
important structures could be lost in the estimate. To tackle this problem, we reformulate the 
reconstruction as a convex optimisation problem and solve it iteratively with an algorithm 
based on Hybrid Steepest Descent (HSD) Yamada (2001). 


We define the multiresolution support M, which is determined by the set of detected 
significant coefficients after hypothesis testing: 


M := {(j,k) if dj[k] is declared significant}. (14) 
We formulate the reconstruction problem as a convex constrained minimization problem: 
Arg min || ®7X||1,s-t. 


X>0, (15) 
{vav € M, (®7X);[k] = (87Y);fk], 


where ® denotes the IUWT synthesis operator. 


This problem is solved with the following iterative scheme: the image is initialised by x) = 
0, and the iteration scheme is, for n = 0 to Nmax — 1: 


Š = P [x + Pua" (Y — xX™)] (16) 
x1) — ØST, [TX] (17) 
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where P, denotes the projection on the positive orthant, Pm denotes the projection on the 
multiresolution support M: 


d;[k] if (j,k) € M, 


Pmdjlk| = { 0 otherwise (18) 
and ST), the soft-thresholding with threshold A: 
= sign(d) (|d| — An) if |d| > Àn, 
STa la] = { 0 otherwise ` (19) 


We chose a decreasing threshold Àn = ei, n= 1,2, $ Nma 


The final estimate of the Poisson intensity is: Â = X(Nmax) | Algorithm 1 summarizes the main 
steps of the MS-VSTS + IUWT denoising algorithm. 


Algorithm 1 MS-VSTS + IUWT Denoising 


Require: data aj := Y, number of iterations Nmax, threshold x 
Detection 

: for j = 1 to J do 

Compute a; and d; using (6). 

Hard threshold |d;[k]| with threshold xo; and update M. 

: end for 
Estimation 

5: Initialize X(°) = 0, Ap = 1. 

6: forn = 0 to Nmax — 1 do 

7 X= P,[xX™ + @PyoeT(y—x™)]. 

8 

9 


eNe 


X("+1) = ST}, [TX]. 
Nmax— 1 
Ant = Nea (nt) ), 
10: end for 
11: Get the estimate Â = X(Nmax), 


4.2 Multi-resolution support adaptation 


When two sources are too close, the less intense source may not be detected because of the 
negative wavelet coefficients of the brightest source. To avoid such a drawback, we may 
update the multi-resolution support at each iteration. The idea is to withdraw the detected 
sources and to make a detection on the remaining residual, so as to detect the sources which 
may have been missed at the first detection. 


At each iteration n, we compute the MS-VSTS of X‘"). We denote a”) [k] the stabilised 


coefficients of X‘"). We make a hard thresholding on (d ilk -— a” [k]) with the same thresholds 


as in the detection step. Significant coefficients are added to the multiresolution support M. 


The main steps of the algorithm are summarized in Algorithm 2. In practice, we use 
Algorithm 2 instead of Algorithm 1 in our experiments. 
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Algorithm 2 MS-VSTS + IUWT Denoising + Multiresolution Support Adaptation 


Require: data aj := Y, number of iterations Nmax, threshold « 
Detection 
1: for j = 1 to J do 
2: Compute aj and dj using (6). 
3: Hard threshold |d;[k]| with threshold xo; and update M. 
4: end for 
Estimation 
: Initialize X) = 0, Ag = 1. 
: for n = 0 to Nmax — 1 do 
X = P [xX + @PyoT(y—x”)]. 
x("+1) = @ST, [TX]. 
Compute the MS-VSTS on X) to get the stabilised coeffcients ae 


10: Hard threshold |d;[k] — d?” [k]| and update M. 


Nmax—(n+1 
11: Anu = Nra (FN 
12: end for 


13: Get the estimate Â = X(Nmax), 


Do o MSI GT 


4.3 MS-VSTS + curvelets 


Insignificant coefficients are zeroed by using the same hypothesis testing framework as in the 
wavelet scale. At each wavelet scale j and ridgelet band k, we make a hard thresholding on 
curvelet coefficients with threshold KO} k, K = 3,4 or 5. Finally, a direct reconstruction can be 
performed by first inverting the local ridgelet transforms and then inverting the MS-VST + 
IUWT (Equation (9)). An iterative reconstruction may also be performed. 


Algorithm 3 summarizes the main steps of the MS-VSTS + Curvelets denoising algorithm. 


Algorithm 3 MS-VSTS + Curvelets Denoising 


1: Apply the MS-VST + IUWT with J scales to get the stabilized wavelet subbands dj. 

2: Set By = Byin. 

3: for j = 1 to J do 

4: Partition the subband d j with blocks of side-length Bj and apply the digital ridgelet 
transform to each block to obtain the stabilized curvelet coefficients. 


5: ifj modulo 2 = 1 then 

6: Bisa = 2B j 

7: else 

8: By = Bj 

9 end if 
10: HTs on the stabilized curvelet coefficients. 
11: end for 


12: Invert the ridgelet transform in each block before inverting the MS-VST + IUWT. 
4.4 Experiments 


The method was tested on simulated Fermi data. The simulated data are the sum of a 
Milky Way diffuse background model and 1000 gamma ray point sources. We based our 
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Galactic diffuse emission model intensity on the model g/l_iem_v02 obtained at the Fermi 
Science Support Center Myers (2009) . This model results from a fit of the LAT photons 
with various gas templates as well as inverse Compton in several energy bands. We used 
a realistic point-spread function for the sources, based on Monte Carlo simulations of the LAT 
and accelerator tests, that scale with energy approximately as 0.8(E/ 1GeV)~°8 degrees (68% 
containment angle). The positions of the 205 brightest sources were taken from the Fermi 
3-month source list Abdo et al. (2009). The positions of the 795 remaining sources follow 
the LAT 1-year Point Source Catalog Myers (2010) source distribution: each simulated source 
was randomly sorted in a box in Galactic coordinates of Al=5° and Ab=1° around a LAT 1-year 
catalog source. We simulated each source assuming a power-law dependence with its spectral 
index given by the 3-month source list and the first year catalog. We used an exposure of 
3.10!9s.cm2 corresponding approximatively to one year of Fermi all-sky survey around 1 GeV. 
The simulated counts map shown in this section correspond to photons energy from 150 MeV 
to 20 GeV. 


Fig. 6 compares the result of denoising with MS-VST + IUWT (Algorithm 1), MS-VST + 
curvelets (Algorithm 3) and Anscombe VST + wavelet shrinkage on a simulated Fermi map. 
Fig. 7 shows the results on one single face of the first scale of the HEALPix pixelization(angular 
extent: 7/3sr). As expected from theory, the Anscombe method produces poor results 
to denoise Fermi data, because the underlying intensity is too weak. Both wavelet and 
curvelet denoising on the sphere perform much better. For this application, wavelets are 
slightly better than curvelets (SNRwavetets = 65-84B, SNReurvelets = 37.3dB, SNR(dB) = 
20 log (Tsignat/ Tnoise)). As this image contains many point sources, this result is expected. 
Indeed wavelets are better than curvelets to represent isotropic objects. 


5. Application to inpainting and source extraction 
5.1 Milky way diffuse background study: denoising and inpainting 


In order to extract from the Fermi photon maps the Galactic diffuse emission, we want to 
remove the point sources from the Fermi image. As our HSD algorithm is very close to the 
MCA (Morphological Component Analysis) algorithm Starck et al. (2004), an idea is to mask 
the most intense sources and to modify our algorithm in order to interpolate through the gaps 
exactly as in the MCA-Inpainting algorithm Abrial et al. (2007). This modified algorithm can 
be called MS-VSTS-Inpainting algorithm. What we want to do is to remove the information 
due to point sources from the maps, in order to keep only the information due to the 
galactic background. The MS-VSTS-Inpainting algorithm interpolates the missing data 
to reconstruct a map of the galactic background, which can now be fitted by a theoretical 
model. The interpolation uses the sparsity of the data in the wavelet domain: the gaps are 
filled so that the result is the sparsiest possible in the wavelet domain. 


The problem can be reformulated as a convex constrained minimization problem: 
Arg min || ®7X||1,s.t. 


x>0, (20) 
ce € M, (TTX); [k] = (®TY) [kK], 


www.intechopen.com 


Poisson Noise Removal in Spherical Multichannel Images: Application to Fermi Data 175 


Simulated Formi dala without noise Simaleted Fermi data with Poisson noise 


Sa a o id C—O 


Simulated Fermi date dennised with Anscombe + envelet shrinknge Simulated Fermi data denoised wilk M3-VST3 + Corvelete 
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oo 
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Fig. 6. Top Left: Fermi simulated map without noise. Top Right: Fermi simulated map with 
Poisson noise. Middle Left: Fermi simulated map denoised with Anscombe VST + wavelet 
shrinkage. Middle Right: Fermi simulated map denoised with MS-VSTS + curvelets 
(Algorithm 3). Bottom Left: Fermi simulated map denoised with MS-VSTS + IUWT 
(Algorithm 1) with threshold 50;. Bottom Right: Fermi simulated map denoised with 
MS-VSTS + IUWT (Algorithm 1) with threshold 3g;. Pictures are in logarithmic scale. 


where IT is a binary mask (1 on valid data and 0 on invalid data). 


The iterative scheme can be adapted to cope with a binary mask, which gives: 


Š = P [x + @Py@ MY —x™)], (21) 
x("+1) = ØST, [Ð]. (22) 


The thresholding strategy has to be adapted. Indeed, for the inpainting task we need to have 
a very large initial threshold in order to have a very smooth image in the beginning and to 
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B W` 

Fig. 7. View of a single HEALPix face (angular extent: 7t /3sr) from the results of Figure 6. Top 
Left: Fermi simulated map without noise. Top Right: Fermi simulated map with Poisson 
noise. Middle Left: Fermi simulated map denoised with Anscombe VST + wavelet shrinkage. 
Middle Right: Fermi simulated map denoised with MS-VSTS + curvelets (Algorithm 3). 
Bottom Left: Fermi simulated map denoised with MS-VSTS + IUWT (Algorithm 1) with 
threshold 50;. Bottom Right: Fermi simulated map denoised with MS-VSTS + IUWT 
(Algorithm 1) with threshold 3g;. Pictures are in logarithmic scale. 
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Strrulated Fermi map with missing ate Simulated Fermi dete denoised and inpainted by MS-VSTS+IU8T 


Fig. 8. MS-VSTS - Inpainting. Left: Fermi simulated map with Poisson noise and the most 
luminous sources masked. Right: Fermi simulated map denoised and inpainted with 
wavelets (Algorithm 4). Pictures are in logarithmic scale. 


refine the details progressively. We chose an exponentially decreasing threshold: 


Nmax=1 


An = Ail Nmax T) — 1), n= 1,2,- , Nmax (23) 


where Amax = max(®"X). 
Algorithm 4 MS-VST + IUWT Denoising + Inpainting 


Require: data ap := Y, mask II, number of iterations Nmax, threshold x. 
Detection 
: for j = 1 to J do 
Compute a; and dj using (6). 
Hard threshold |d;[k]| with threshold xo; and update M. 
: end for 
Estimation 
5: Initialize X(0) = 0, Ag = Amax- 
6: forn = 0 to Nmax — 1 do 
7. 
8 


eNe 


X = P [x + #PuSTT(Y — X™)]. 
x+) = ost [7X]. 


(Mme) 
9: Anti = Amax (2 Nae 1) 
10: end for 


11: Get the estimate A = X(Nmax), 


Experiment 


We applied this method on simulated Fermi data where we masked the 500 most luminous 
sources. (with the highest photon per pixel flux) The other sources are not intense enough 
to be differencied from the background. 


The results are on Figure 8. The MS-VST + IUWT + Inpainting method (Algorithm 4) 
interpolates the missing data very well. Indeed, the missing part can not be seen anymore 
in the inpainted map, which shows that the diffuse emission component has been correctly 
reconstructed. 
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5.2 Source detection: denoising and background modeling 
5.2.1 Method 


In the case of Fermi data, the diffuse gamma-ray emission from the Milky Way, due to 
interaction between cosmic rays and interstellar gas and radiation, makes a relatively intense 
background. We have to extract this background in order to detect point sources. This diffuse 
interstellar emission can be modelled by a linear combination of gas templates and inverse 
compton map. We can use such a background model and incorporate a background removal 
in our denoising algorithm. 


We denote Y the data, B the background we want to remove, and a” [k] the MS-VSTS 
coefficients of B at scale j and position k. We determine the multi-resolution support by 
comparing |d;[k] — a [k]| with Koj. 


We formulate the reconstruction problem as a convex constrained minimization problem: 
Arg min || ®7X||1,s.0. 


X>0, (24) 
{van E€ M, (®TX) [kK] = (87(Y — B)) [A], 


Then, the reconstruction algorithm scheme becomes: 


X = P [x + Puo (y—B-x”)], (25) 
x("+)) = @ST, [TX]. (26) 


The algorithm is illustrated by the theoretical study in Figure 9. We denoise Poisson data 
while separating a single source, which is a Gaussian of standard deviation equal to 0.01, 
from a background, which is a sum of two Gaussians of standard deviation equal to 0.1 and 
0.01 respectively. 


Algorithm 5 MS-VSTS + IUWT Denoising + Background extraction 


Require: data ap := Y, background B, number of iterations Nmax, threshold x. 
Detection 
1: for j = 1 to J do 
2: Compute aj and dj using (6). 
3: Hard threshold (d;[k] — a [k]) with threshold xo; and update M. 
4: end for 
Estimation 
5: Initialize X®) = 0, Ag = 1. 
6: forn = 0 to Nmax — 1 do 
77 X= P [X +4 @P\yoT(y-B-x™)]. 
8: X("1) = OST, [TX]. 
9 Nmax—("+1)_ 


max 


Anti al 
10: end for 
11: Get the estimate A = X(Nmax), 
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Fig. 9. Theoretical testing for MS-VSTS + IUWT denoising + background removal algorithm 
(Algorithm 5). View on a single HEALPix face. Top Left: Simulated background : sum of two 
Gaussians of standard deviation equal to 0.1 and 0.01 respectively. Top Right: Simulated 
source: Gaussian of standard deviation equal to 0.01. Bottom Left: Simulated poisson data. 
Bottom Right: Image denoised with MS-VSTS + IUWT and background removal. 


Like Algorithm 1, Algorithm 5 can be adapted to make multiresolution support adaptation. 


5.2.2 Experiment 


We applied Algorithms 5 on simulated Fermi data. To test the efficiency of our method, 
we detect the sources with the SExtractor routine Bertin & Arnouts (1996), and compare the 
detected sources with the input source list to get the number of true and false detections. 
Results are shown on Figures 10 and 11. The SExtractor method was applied on the first 
wavelet scale of the reconstructed map, with a detection threshold equal to 1. It has been 
chosen to optimise the number of true detections. SExtractor makes 593 true detections and 
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Fig. 10. Top Left: Simulated background model. Top Right: Simulated Gamma Ray sources. 
Middle Left: Simulated Fermi data with Poisson noise. Middle Right: Reconstructed Gamma 
Ray Sources with MS-VSTS + IUWT + background removal (Algorithm 5) with threshold 50;. 
Bottom: Reconstructed Gamma Ray Sources with MS-VSTS + IUWT + background removal 
(Algorithm 5) with threshold 3g;. Pictures are in logarithmic scale. 


71 false detections on the Fermi simulated map restored with Algorithm 2 among the 1000 
sources of the simulation. On noisy data, many fluctuations due to Poisson noise are detected 
as sources by SExtractor, which leads to a big number of false detections (more than 2000 in 
the case of Fermi data). 


Sensitivity to model errors 


As it is difficult to model the background precisely, it is important to study the sensitivity of 
the method to model errors. We add a stationary Gaussian noise to the background model, we 
compute the MS-VSTS + IUWT with threshold 3g; on the simulated Fermi Poisson data with 
extraction of the noisy background, and we study the percent of true and false detections with 
respect to the total number of sources of the simulation and the signal-noise ratio (SNR(dB) = 
20 log (signal/Cnoise)) versus the standard deviation of the Gaussian perturbation. Table 3 


www.intechopen.com 


Poisson Noise Removal in Spherical Multichannel Images: Application to Fermi Data 181 


Fig. 11. View of a single HEALPix face (angular extent: 7t /3sr) from the results of 

Figure 10.Top Left: Simulated background model. Top Right: Simulated Gamma Ray sources. 
Middle Left: Simulated Fermi data with Poisson noise. Middle Right: Reconstructed Gamma 
Ray Sources with MS-VSTS + IUWT + background removal (Algorithm 5) with threshold 507. 
Bottom: Reconstructed Gamma Ray Sources with MS-VSTS + IUWT + background removal 
(Algorithm 5) with threshold 3g;. Pictures are in logarithmic scale. 


www.intechopen.com 


182 Wavelet Transforms and Their Recent Applications in Biology and Geoscience 


Model error std dev|% of true detect|% of false detect 
0 59.3% 7.1% 23.8 
10 57.0% 11.0% 23.2: 
20 53.2% 18.9% 22.6 
30 49.1% 43.5% 21.7 
40 42.3% 44.3% 21.0 
50 34.9% 39.0% 20.3 
60 30.3% 37.5% 19.5 
70 25.0% 34.6% 18.9 
80 23.0% 28.5% 18.7 
90 23.6% 27.1% 18.3 


Table 3. Percent of true and false detection and signal-noise ratio versus the standard 
deviation of the Gaussian noise on the background model. 


shows that, when the standard deviation of the noise on the background model becomes of 
the same range as the mean of the Poisson intensity distribution (Amean = 68.764), the number 
of false detections increases, the number of true detections decreases and the signal noise ratio 
decreases. While the perturbation is not too strong (standard deviation < 10), the effect of the 
model error remains low. 


6. Extension to multichannel data 
6.1 Gaussian noise 


6.1.1 2D-1D Wavelet Transform on the sphere 


We propose a denoising method for 2D - 1D data on the sphere, where the two first dimensions 
are spatial (longitude and latitude) and the third dimension is either the time or the energy. 
We need to analyze the data with a non-isotropic wavelet, where the time or energy scale is 
not connected to the spatial scale. An ideal wavelet function would be defined by: 


(8, ,t) = peP 0, pp E) (27) 


where 9) is the spatial wavelet and yl) is the temporal (or energy) wavelet. In the 
following, we will consider only isotropic and dyadic spatial scales, and we denote jı the 
spatial resolution index (i.e. scale 2/!), jọ the time (or prey Seed index. We thus 


define j scaled spatial and temporal (or energy) wavelets p ) (o, 9) = J 1 (8,9) ( £ i of ) 
and pt") = 5 zp). 


Hence, we derive the wavelet coefficients Wj ja [ko, ko, ky] from a given data set D (kg and ko 
are spatial index and kz a time (or energy) index. In continuous coordinates, this amounts to 
the formula 


Wj j lko, ko, kt] = GR Hye D(0, g,t) 


0 — kg P- ko (t) kt = z(t) 
J OR y q )dxdydz = Dx ple = pt (0, @, t) 


(28) 


x (FP) ( 
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where x is the convolution and #(t) = (—t). 


6.1.2 Fast undecimated 2D-1D decomposition/reconstruction 


In order to have a fast algorithm for discrete data, we use wavelet functions associated to filter 
banks. Hence, our wavelet decomposition consists in applying first a IUWT on the sphere for 
each frame kz. Using the spherical IUWT, we have the reconstruction formula: 


J1 
D{ko, ko, ki] =a), (ko, kọ] + L Wi ko, ko, kt], Vki (29) 
n=l 


where J; is the number of spatial scales. To have simpler notations, we replace the two spatial 
indexes by a single index k, which corresponds to the pixel index: 
yl 
D{ky, kt] = ay, [ky] + L Wj [kr, kt], Vki (30) 
j=1 


Then, for each spatial location k, and for each 2D wavelet scale jı, we apply a 1D wavelet 
transform along t on the spatial wavelet coefficients Wj, [kr, ke] such that 


h 


wj [kr kt] = wi, Jy [kr kt] + D Wi, ja [Krs kt], Y (kr, ke) (31) 
p= 


where j2 is the number of scales along t.The same processing is also applied on the coarse 
spatial scale aj, [ky,k;] and we have 


In 
ay, [ky, ki] =A), Jo [ky, ki] + 3 Why, jo [ky, ke], V(ky, kt) (32) 
p= 
Hence, we have a 2D-1D spherical undecimated wavelet representation of the input data D: 
J Jo Jy——J2 
Dik, kt] = aj, Jy [ky, ke] + B Wih [ky, ke] T D Whja [ky, ke] + L y Wija [ky, ke] (33) 
h=1 p= h=1p=l 

From this expression, we distinguish four kinds of coefficients: 
e Detail-Detail coefficients (j4 < J4 and jz < J2): 


w; jplke, kt] = (6 — mp) x (AYR? xap ake, -] — YB? xaj [ke ]) (34) 


e Approximation-Detail coefficients (j4 = J, and jz < J2): 


wp jalkr, ke] = WY”) xap [kr] — WY) xap [kr] (35) 
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e Detail-Approximation coefficients (j4 < J, and j2 = J2): 
wjp [krki] = hB xap —a lke] = AB * aj lhe (36) 
e Approximation-Approximation coefficients (jı = J, and j2 = Jo): 
ajj lkr ke] = iB xap [hry 67) 


6.1.3 Multichannel Gaussian denoising 


As the spherical 2D-1D undecimated wavelet transform just described is fully linear, a 
Gaussian noise remains Gaussian after transformation. Therefore, all thresholding strategies 
which have been developed for wavelet Gaussian denoising are still valid with the spherical 
2D-1D wavelet transform. Denoting TH the thresholding operator, the denoised cube in the 
case of additive white Gaussian noise is obtained by: 


Djk, ky] = = Aj, Jy [ky, ke] ]+ 2 TH(w Wih [k,, k] ) 
j=1 
Jo h b 
+ 2 TH(wy, jo[kr, ke]) + D> } TH(w;,j [kr ke]) (38) 
p= j=1j:2=1 
A typical choice of TH is the hard thresholding operator, i.e. 
Soif 
TH(x) = T Fidar (39) 


The threshold qt is generally chosen between 3 and 5 times the noise standard deviation. 


6.2 Poisson Noise 
6.2.1 Multi-scale variance stabilzing transform 


To perform a Poisson denoising, we have to plug the MS-VST into the spherical 2D-1D 
undecimated wavelet transform. Again, we distinguish four kinds of coefficients that take 
the following forms: 


e Detail-Detail coefficients (j1 < J; and jz < J2): 
wi lkr ki] = (8 — ħin) * (Taha gta ABP xa ilk Th— Trea [2 PP xa [hes ]]) (40 
jitt Rt] = 11D x ( ALL pti 1D * aj, 11 vr J] in iI ip * aj, | t7 I) ( ) 
e Approximation-Detail coefficients (jı = J, and jz < Jo): 


=] 


wp jplkr Rt] = Th ja ABP xap fkr 1] — Tp jn A? x a, [kr I] (41) 


e Detail-Approximation coefficients (j1 < Jı and j2 = Jo): 


Wj polkre ke] = Taip AP) « aj, 1 lhe, J] — Thy pA xaj [kr] (42) 
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e Approximation-Approximation coefficients (jj = J, and j2 = J): 
Aj Jalker ke] = hi xap [hry (43) 


Hence, all 2D-1D wavelet coefficients Wiji, j, are NOW stabilized, and the noise on all these 
wavelet coefficients is Gaussian with known scale-dependent variance that depends solely 
on h. Denoising is however not straightforward because there is no explicit reconstruction 
formula available because of the form of the stabilization equations above. Formally, the 
stabilizing operators Tjj and the convolution operators along the spatial and temporal 
dimensions do not commute, even though the filter bank satisfies the exact reconstruction 
formula. To circumvent this difficulty, we propose to solve this reconstruction problem by 
using an iterative reconstruction scheme. 


6.2.2 Detection-reconstruction 


As the noise on the stabilized coefficients is Gaussian, and without loss of generality, we let its 
standard deviation equal to 1, we consider that a wavelet coefficient Wija [k;, kt] is significant, 
i.e., not due to noise, if its absolute value is larger than a critical threshold t, where T is 


typically between 3 and 5. 


The multiresolution support will be obtained by detecting at each scale the significant 
coefficients. The multiresolution support for j; < J; and j2 < J2 is defined as: 


1 if w; ;, [kr, kilis significant 


as = hide 
Misji [kr ke] e otherwise ma 


We denote W the spherical 2D-1D undecimated wavelet transform described above, and R 
the inverse wavelet transform. We want our solution X to preserve the significant structures of 
the original data by reproducing exactly the same coefficients as the wavelet coefficients of the 
input data Y, but only at scales and positions where significant signal has been detected. At 
other scales and positions, we want the smoothest solution with the lowest budget in terms 
of wavelet coefficients.Furthermore, as Poisson intensity functions are positive by nature, a 
positivity constraint is imposed on the solution. It is clear that there are many solutions 
satisfying the positivity and multiresolution support consistency requirements, e.g. Y itself. 
Thus, our reconstruction problem based solely on these constraints is an ill-posed inverse 
problem that must be regularized. Typically, the solution in which we are interested must be 
sparse by involving the lowest budget of wavelet coefficients. Therefore our reconstruction is 
formulated as a constrained sparsity-promoting minimization problem that can be written as 
follows 


f 5 MWX = MWY 
min || WVX||_ subject to { X>0 (45) 
where ||- || is the Lj-norm playing the role of regularization and is well known to promote 


sparsityDonoho (2004). This problem can be solved efficently using the hybrid steepest 
descent algorithm Yamada (2001)Zhang et al. (2008), and requires about 10 iterations in 
practice. Transposed into our context, its main steps can be summarized as follows: 
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Require: Input noisy data Y, a low-pass filter h, multiresolution support M from the 
detection step, number of iterations Nmax 
: Initialize XO) = MWY = Mwy, 
: forn = 1 to Nmax do 
d= Mwy + (1—-M)wx("-)), 
x") = P,(RSTz, [d]), 
Update the step By = (Nmax — 1)/(Nmax — 1) 
end for 


ON te Ore 


where P, is the projector onto the positive orthant, i.e. P}(x) = max(x,0), STg, is the 
soft-thresholding operator with threshold By, ie. ST, [x] = x — Bnsign(x) if |x| > Bn, and 
0 otherwise. 


The final spherical MSVST 2D-1D wavelet denoising algorithm is the following: 


Require: Input noisy data Y, a low-pass filter h, threshold level t 
1: Spherical 2D-1D MSVST: Apply the spherical 2D-1D-MSVST to the data using (40)-(43). 
2: Detection: Detect the significant wavelet coefficients that are above t, and compute the 
multiresolution support M. 
3: Reconstruction: Reconstruct the denoised data using the algorithm above. 


7. Conclusion 


This chapter presented new methods for restoration of spherical data with noise following a 
Poisson distribution. A denoising method was proposed, which used a variance stabilization 
method and multiscale transforms on the sphere. Experiments have shown it is very efficient 
for Fermi data denoising. Two spherical multiscale transforms, the wavelet and the curvelets, 
were used. Then, we have proposed an extension of the denoising method in order to take 
into account missing data, and we have shown that this inpainting method could be a useful 
tool to estimate the diffuse emission. Then, we have introduced a new denoising method 
on the sphere which takes into account a background model. The simulated data have 
shown that it is relatively robust to errors in the model, and can therefore be used for Fermi 
diffuse background modeling and source detection. Finally, we introduced an extension for 
multichannel data. 
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